//
//---------------------  1D  ---------------------
//
//
// Upwind is not defined for 1D
//

//
//---------------------  2D  ---------------------
//
//       Staggered Mesh for u-vel and v-vel
//
//       0       1       2       3       4          
//
//5      >       >       >       >       > 
//       |       |       |       |       |              
//   ^---+---^---+---^---+---^---+---^---+---^  4       Mesh for scalar fields
//       |       |       |       |       | 
//4      >   o   >   o   >   o   >   o   >               5  x-x-+-x-+-x-+-x-x
//       |       |       |       |       |               4  x o | o | o | o x
//   ^---+---^---+---^---+---^---+---^---+---^  3           +---+---+---+---+
//       |       |       |       |       |               3  x o | o | o | o x
//3      >   o   >   o   >   o   >   o   >                  +---+---+---+---+
//       |       |       |       |       |               2  x o | o | o | o x
//   ^---+---^---+---^---+---^---+---^---+---^  2           +---+---+---+---+
//       |       |       |       |       |               1  x o | o | o | o x
//2      >   o   >   o   >   o   >   o   >               0  x-x-+-x-+-x-+-x-x
//       |       |       |       |       |                  
//   ^---+---^---+---^---+---^---+---^---+---^  1           0 1   2   3   4 5 
//       |       |       |       |       |                 
//1      >   o   >   o   >   o   >   o   >                 o central node
//       |       |       |       |       |                 x boundary node
//   ^---+---^---+---^---+---^---+---^---+---^  0          > u velocity 
//       |       |       |       |       |                 ^ v velocity
//0      >       >       >       >       >    
//   0       1       2       3       4       5       
//                                                  
//                      v(i,j) 
//                  |     n     |               
//                --|-----^-----|--             
//                  |           |               
//                  |           |               
//     u(i-1,j) = w >     o     > e = u(i,j)
//                  |   (i,j)   |
//                  |           |
//                --|-----^-----|--
//                  |     s     |
//                     v(i,j-1)
//               
//
//---------------------  3D  ---------------------
//
template<class T_number, int Dim>
inline 
bool UpwindE<T_number, Dim>::calcCoefficients(const ScalarField &nut) { 
    T_number dyz = dy * dz, dxz = dx * dz, dxy = dx * dy;
    T_number dyz_dx = dyz / dx, dxz_dy = dxz / dy, dxy_dz = dxy / dz;
    T_number ce, cw, cn, cs, cf, cb;
    T_number nutinter;
    T_number dxyz_dt = dx * dy * dz / dt;

    T_number Prt = PrandtlTurbulent;

    for (int i =  bi; i <= ei; ++i)
	for (int j = bj; j <= ej; ++j)
	    for (int k = bk; k <= ek; ++k)
	    {
		ce = u(i  , j, k) * dyz;
		cw = u(i-1, j, k) * dyz;		
		if ( ce > 0 ) ce = 0.0; 
		else          ce = -ce;
		if ( cw <= 0 ) cw = 0.0;

		cn = v(i, j  , k) * dxz;
		cs = v(i, j-1, k) * dxz;		
		if ( cn > 0 ) cn = 0.0;
		else          cn = -cn;
		if ( cs <= 0 ) cs = 0.0; 

		cf = w(i, j, k) * dxy;
		cb = w(i, j, k-1) * dxy;		
		if ( cf > 0 ) cf = 0.0;
		else          cf = -cf;
		if ( cb <= 0 ) cb = 0.0; 

/*********************************************************************
// Actually, nut is calculated on center of volumes!!!
//		
// kappa_turbulent = nu_turbulent / Pr_turbulent.
// Because nut is calculated on nodes, we have to interpolate to the
// center of volume.
//
		nutinter = 0.125 * ( nut(i  ,j  ,k  ) + nut(i  ,j  ,k-1) +
				     nut(i  ,j-1,k-1) + nut(i  ,j-1,k  ) +
				     nut(i-1,j  ,k  ) + nut(i-1,j  ,k-1) +
				     nut(i-1,j-1,k-1) + nut(i-1,j-1,k  ) );

		aE (i,j,k) = (1 + nutinter / Prt) * dyz_dx + ce;
		aW (i,j,k) = (1 + nutinter / Prt) * dyz_dx + cw;
		aN (i,j,k) = (1 + nutinter / Prt) * dxz_dy + cn;
		aS (i,j,k) = (1 + nutinter / Prt) * dxz_dy + cs;
		aF (i,j,k) = (1 + nutinter / Prt) * dxy_dz + cf;
		aB (i,j,k) = (1 + nutinter / Prt) * dxy_dz + cb;
*********************************************************************/

		aE (i,j,k) = (1 + nut(i,j,k) / Prt) * dyz_dx + ce;
		aW (i,j,k) = (1 + nut(i,j,k) / Prt) * dyz_dx + cw;
		aN (i,j,k) = (1 + nut(i,j,k) / Prt) * dxz_dy + cn;
		aS (i,j,k) = (1 + nut(i,j,k) / Prt) * dxz_dy + cs;
		aF (i,j,k) = (1 + nut(i,j,k) / Prt) * dxy_dz + cf;
		aB (i,j,k) = (1 + nut(i,j,k) / Prt) * dxy_dz + cb;
		aP (i,j,k) = aE (i,j,k) + aW (i,j,k) + 
		             aN (i,j,k) + aS (i,j,k) +
		             aF (i,j,k) + aB (i,j,k) + dxyz_dt;
//		+ (ce - cw) + (cn - cs) + (cf - cb);

		sp (i,j,k) = T(i,j,k) * dxyz_dt ;
	    }
//    cout << "\n Dim :" << dimension << " function : a_EW_3D ";
    applyBoundaryConditions();

    return 1;
}
















